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METHOD AND APPARATUS FOR PRODUCING A REPRESENTATION OF 
A MEASURABLE PROPERTY WHICH VARIES IN TIME AND SPACE, FOR 
PRODUCING AN IMAGE REPRESENTING CHANGES IN RADIOACTIVITY 
IN AN OBJECT AND FOR ANALYZING TOMOGRAPHY SCAN IMAGES 

5 

RELATED APPLICATIONS 

This application claims priority from United States provisional application 
serial No. 60/107,335, filed Novembers, 1998. 



1 0 FIELD OF THE INVENTION 

The present invention relates to dynamic data analysis, and more particularly 
to methods and apparatus for producing a representation of a measurable 
property which varies in time and space, for producing an image representing 
changes in radioactivity in an object and for analyzing tomography scan 
15 images. 



BACKGROUND OF THE INVENTION 

In recent years, nuclear medicine procedures have been increasingly used to 
trace distributions of various substances in the human body. Typically, 

20 radioisotope-tagged pharmaceuticals (radiopharmaceuticals) are administered 

to a subject, which may be a human or an animal, for example. Common 
imaging techniques, such as Single Photon Emission Computed Tomography 
(SPECT), for example, are then used to measure radiation from the 
radiopharmaceuticals to produce data and/or images representing the three- 

25 dimensional distribution of the radiopharmaceuticals in the subject's body or in 



a particular organ or region thereof. Scans of this nature may be used for 
clinical diagnosis of sick patients, for example. 

Current SPECT visualization is generally perfomied under the assumption 
that the spatial distribution of the radiopharmaceuticals in the subject's body is 
static throughout the duration of the scan, typically 20-25 minutes. However, 
physiological processes in the human body are usually dynamic, therefore 
many organs, such as the kidney and heart, for example, tend to show 
significant changes in activity distribution over time due to uptake (Vashin") 
or due to "washout" of the radiopharmaceuticals. 

Moreover, it is presently thought that the rates of infusion and subsequent 
extraction of such radiopharmaceuticals from the tissue provide good 
indications of organ function. Therefore, while current SPECT technology and 
methods are useful for measuring a static distribution of a substance, they are 
not suitable for measuring dynamic physiological processes, such as kidney 
or heart function, for example, which require dynamic data representing the 
"washin" or "washout" rates of the substance. 

In other words, measurement of physiological processes requires four- 
dimensional dynamic data representing both the spatial and temporal 
distributions of the radiophannaceutical. 

Other technologies exist for producing dynamic data of this general nature, 
however, such technologies suffer from a number of disadvantages. 

One technique involves fast rotations around the subject with a SPECT 
camera. However, due to the relatively high speed at which the camera 
moves, very few counts are received at any given camera location. This 
method therefore produces statistically unreliable data, resulting in poor 
spatial resolution. Additionally, this technique is generally incompatible with 
the majority of transmission scan-based attenuation correction methods. 
Equipment for perfonning this technique is rarely available in hospitals. 



Ring SPECT devices may also suffer from statistically unreliable data. 
Typically, each of the detectors forming the ring receives a relatively low 
number of counts, resulting in a low signal to noise ratio. Ring SPECT 
devices are also expensive and are typically not available in hospitals, being 
5 used mainly for research. 

Similarly, Positron Emission Tomography (PET) devices are also prohibitively 
expensive, and are not available in many hospitals. PET equipment is much 
more complex than SPECT equipment, and the PET technique requires a 
cyclotron in proximity to or in the hospital in which the PET equipment is 
10 housed. In addition, PET and SPECT use different isotopes and 

radiopharmaceuticals, and are therefore complementary. In addition, data 
analysis may be cumbersome with the PET technique, wherein typically 64 or 
128 projections or input images are required to produce each reconstructed 
image. 

15 Planar imaging is commonly available in hospitals. However, this technique 

generally involves acquiring a series of fast images, typically 1-10 seconds, 
which sometimes results in a low signal-to-noise ratio. The camera generally 
remains at a fixed position, resulting in limited spatial resolution, an inability to 
produce three-dimensional images of the object of interest and an inability to 

20 perform proper attenuation correction. 

More generally, existing dynamic analysis methods generally assume that the 
activity of each three-dimensional pixel or "voxel" of the area of interest 
behaves according to a particular functional model, such as exponential 
decay or biexponential decay, for example. One such method, Factor 

25 Analysis of Dynamic Structures (FADS), proposes a "menu" of possible 

dynamic profiles or basis functions, sometimes referred to as "factors", which 
may include exponential or biexponential functions for example. FADS 
attempts to fit a linear combination of such basis functions to the measured 
data. However, the exponential or other functional models imposed on the 

30 data may not apply in some situations, and in such a case this technique may 



fail. This method may also fail due to inherent numerical instabilities if the 
number of dynamic pixels is high. In addition, existing dynamic analysis 
methods are cumbersome and time-consuming, and often require as much as 
one to two hours or more of computing time to complete the data analysis 
using contemporary desktop computers. 

Accordingly, there is a need for a way to produce a representation of a 
measurable property which varies in time and space. More particularly, there 
is a need for ways to produce an image representing changes in radioactivity 
in an object and to use standard nuclear medicine equipment currently 
available in most hospitals, to quickly produce dynamic data representing 
physiological processes over time, without sacrificing the resolution or 
reliability of the data, and without requiring the data to satisfy any pre- 
determined form of function. 



SUMMARY OF THE INVENTION 

The present invention addresses the above need by providing a computer 
implemented method of producing a representation of a measurable property 
which varies in time and space. A plurality of sets of values representing 
measurements of the property across an object are received, each set being 
associated with a respective measurement time. A plurality of sets of values 
representing the property at a plurality of locations throughout the object at 
the respective measurement times are then produced. This is achieved by 
minimizing a figure of merit function which relates the values representing the 
measurements with the values representing the property at the plurality of 
locations, subject to a shape constraint imposed on the values representing 
the property at the plurality of locations throughout the object. 

Particular embodiments of the present invention provide a novel and inventive 
method of reconstruction using Shape Constrained Least Squares (SCLS) 
analysis. The figure of merit function is minimized without making any 



assumptions as to the precise functional form of tlie time variation of the 
measurable property, but rather, assuming only a shape constraint. 

One particularly advantageous embodiment of the invention provides a 
Dynamic Single Photon Emission Computed Tomography (DSPECT) method 
and apparatus. Conventional SPECT imaging devices, such as those found 
in many modern hospitals, may be used to produce data representing 
tomography scan images of the object, which may be a functioning heart or 
kidney contained within a living person, for example, to which a 
radiopharmaceutical has been administered. In such a case, SCLS analysis 
according to the present invention produces dynamic data representing the 
time variation of the three-dimensional radioactivity distribution across the 
object, throughout the duration of the tomography scan, typically 20-25 
minutes. This dynamic information may provide indications of organ function, 
for example, which cannot be obtained from conventional SPECT techniques 
in which the data is assumed to be static. 

In this embodiment, the present invention therefore permits use of 
conventional equipment commonly available in hospitals to produce dynamic 
information previously obtainable only by use of expensive and/or rare 
equipment such as PET or Ring SPECT devices, for example. In addition, the 
use of conventional SPECT equipment results in a relatively high signal to 
noise ratio compared to certain types of equipment previously used to 
produce such dynamic information. 

In contrast with analysis techniques previously used to produce dynamic data, 
such as Non-linear Least Squares (NLS) analysis, for example, the inventive 
SCLS analysis method disclosed herein provides for relatively fast data 
analysis, producing the dynamic data representing the time variation of the 
physical property in as little as ten minutes, for example, as compared to 
typical computing times of one to two hours or more with NLS analysis 
techniques. More importantly, the present invention does not assume any 
precise functional form of the dynamic data, such as the exponential or 
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biexponential models commonly used in NLS analysis, which may not 
necessarily apply in all situations. 

Examples of clinical applications of this embodiment of the invention may 
include renal studies using the dynamic tracers Tc-99mMAG3 and Tc- 
5 99mDTPA used in planar renal imaging, Tc-99m labeled teboroxime studies 

used to investigate myocardial blood flow for detection of ischemia or 
infarction, myocardial viability studies using 1-123 labeled fatty acids, or brain 
blood-flow imaging using Xe-133, Xe-127 or 1-123, for example. The 
biological half-lives of any of these dynamic tracers is within the time scale of 
1 0 a typical scan of 20-25 minutes. 

An apparatus for performing the foregoing method may include a receiver for 
receiving the plurality of sets of values representing the measurements of the 
property, and a processor circuit in communication with the receiver. The 
processor circuit is programmed to produce the plurality of sets of values 
15 representing the property at the plurality of locations throughout the object, 

according to the methods indicated above and herein. 

In accordance with another aspect of the invention, there is provided a 
computer readable medium for providing computer readable instructions for 
directing a programmable device to implement the method indicated above. 
20 Similarly, another aspect of the invention involves providing a computer data 

signal embodied in a carrier wave, the data signal including code segments 
for directing a programmable device to implement the methods indicated 
above and herein. 

According to yet another aspect of the invention, there is provided a method 
25 of producing an image representing changes of radioactivity in an object. The 

method includes receiving a plurality of sets of values representing 
tomography scan images across the object at respective measurement times, 
each set being associated with a respective measurement time. The method 
further includes producing a plurality of sets of values representing 
30 radioactivity at a plurality of locations throughout the object at the respective 
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measurement times, by minimizing a figure of merit function relating the 
values representing the tomography scan images with the values representing 
radioactivity, with a shape constraint imposed on the values representing 
radioactivity. The method also includes producing a visual representation of 
5 the object in response to the plurality of sets of values representing 

radioactivity, the visual representation including a representation of 
radioactivity over time. An apparatus for performing the method includes a 
receiver for receiving the plurality of sets of values representing the 
tomography scan images, and a processor circuit programmed to produce the 

10 plurality of sets of values representing radioactivity, as indicated. The 

processor circuit is also programmed to produce a time varying graphical 
representation of the object in response to the plurality of sets of values 
representing radioactivity, and the apparatus also Includes a display 
responsive to the graphical representation for producing a visual 

15 representation of the object in response to the graphical representation. 

According to an additional aspect of the invention there is provided a method 
of using a processor to analyze data signals representing tomography scan 
images of an organic object. The method includes receiving data 
representing successive tomography scan images of the object, and 

20 performing calculations by imposing an inequality constraint to determine 

dynamic data values from the data. Each of the dynamic data values 
represents a physical property of the object at a respective corresponding one 
of a plurality of voxels of the object at a respective corresponding time. The 
method also includes producing a representation of the dynamic data values, 

25 representing the physical property at the voxels at the times. 

Other aspects and features of the present invention will become apparent to 
those ordinarily skilled in the art upon review of the following description of 
specific embodiments of the invention in conjunction with the accompanying 
figures. 
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BRIEF DESCRIPTION OF THE DRAWINGS 



In drawings which illustrate embodiments of the invention, 

Figure 1 is a schematic representation of a system for producing a 
representation of a measurable property which varies in time and 
space, according to a first embodiment of the invention. 

Figure 2 is a fragmented schematic representation of an imaging device 
and an object shown in Figure 1. 

Figure 3 is a block diagram of a programmable device of the system shown 
in Figure 1 . 

Figure 4 is a flow chart of a receive data routine executed by the 
programmable device shown in Figure 3. 

Figure 5 is a flow chart of a dynamic solution routine executed by the 
programmable device shown in Figure 3. 

Figure 6 is a schematic representation of a dynamic data matrix produced 
by the programmable device shown in Figure 3. 



DETAILED DESCRIPTION 

As shown in Figure 1, a system for producing a representation of a measurable 
property which varies in time and space is designated generally by the reference 
character 20. The system includes a data gathering device shown generally at 
22 for producing signals or data values representing the measurable property of 
an object of interest 24 at respective measurement times, and a programmable 
device shown generally at 26. 

In this embodiment, the system 20 includes a system for Dynamic Single Photon 
Emission Computed Tomography (DSPECT). The object 24 includes an 



organic object 28, such as a kidney or a heart, for example. The organic object 
28 is preferably contained within a living organism (not shown) such as a human 
or an animal, for example. 

Also in this embodiment, the programmable device 26 includes a desktop 
computer 34 having a user input device 33, which in this embodiment includes a 
keyboard and a mouse. The computer further has an output device, which in 
this embodiment includes a display device 35 for displaying a visual 
representation of the measurable property. 

Data gathering device 

Generally, the data gathering device 22 produces a plurality of sets of values 
representing measurements of the measurable property across the object 28 at 
respective measurement times. Each set of values is associated with a 
respective one of the measurement times. 

As shown in Figure 1 , in this embodiment the data gathering device 22 includes 
a tomographic imaging device 30 and a control device 32 for moving the 
tomographic imaging device 30 relative to the organic object 28. 

The control device 32 controls the tomographic imaging device 30 to move the 
tomographic imaging device along a 180" circumferential arc 31 about the object 
28 at a radial distance of approximately 25 cm from the centre of the object. 
More particularly, in this embodiment the control device moves the tomographic 
imaging device among ninety discrete positions or stops along the 180' 
circumferential arc 31, with each stop lasting approximately 13 seconds, for a 
total data acquisition time over the complete arc 31 of approximately 20 minutes. 
At each stop, the tomographic imaging device produces values representing a 
tomography scan image of the object. These values are received and stored by 
the programmable device 26, following which the control device moves the 
tomographic imaging device to the next stop, and so on, until the 180° scan has 
been completed. 
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Alternatively, it will be appreciated that the number of stops, duration of each 
stop, arc length or angle of rotation and radial distance may be varied. 
Generally, distances as close as possible to the object may be preferred. The 
time duration of each stop may be tailored to a particular clinical problem, and 
will be related to the temporal resolution of the study. Where the object is an 
organ in a human, it has been found that between 60 and 120 stops along a 
180" arc at a radial distance of 20 to 30 cm from the object is preferable. 

Tomographic imaging device 

As shown in Figure 2, the tomographic imaging device is designated generally 
by the reference character 30. Generally, the tomographic imaging device is 
used to produce successive sets of values representing successive tomography 
scan images across the object 28 at respective measurement times, each set of 
values being associated with a respective measurement time. In this 
embodiment, the tomographic imaging device includes a collimator 36 and a 
detector shown generally at 38. The tomographic imaging device may be 
conventional, however, any tomographic imaging device capable of performing 
the functions described herein may be substituted therefor. In this embodiment, 
the tomographic imaging device includes a Single Photon Emission Computed 
Tomography (SPECT) imaging device, which acts as a tomography scanner for 
measuring radioactivity across the object. Alternatively, the imaging device may 
include a dual or triple head camera to provide better accuracy, or a coincidence 
SPECT camera, for example. Preferably, the tomographic imaging device 30 is 
capable of perfomning a transmission scan to allow for better attenuation 
correction. 

For analytical purposes, the object 28 is defined as comprising a plurality of 
cubical regions called voxels, some of which are shown generally at 40 in Figure 
2. Typically, the object may comprise 30 x 30 x 30 voxels of 1 cm^ each, for 
example, however, for illustrative purposes the size of each voxel is 
exaggerated in Figure 2. 



-11- 

Typically, a radioisotope-labelled pharmaceutical or "radiopharmaceutical" is 
injected or otherwise administered into the object 28. In this embodiment, the 
radioisotope is one that emits gamma radiation, such as Tc-99m, 1-123 or Xe- 
127, for example. More particularly, in this embodiment the radiopharmaceutical 
5 includes Tc-99m labeled teboroxime. However, it will be appreciated that other 

radiopharmaceuticals may be substituted therefor. 

Accordingly, in this embodiment the detector 38 includes a gamma ray detector, 
which in turn includes an array of receptor bins shown generally at 42 in Figure 
2. The detector 38 is approximately 30 x 40 cm, and each of the bins 42 is 
10 approximately 6x6 mm. In this embodiment, therefore, the array comprises 

approximately 50 x 66 receptor bins 42, for a total of approximately 3300 
receptor bins. For illustrative purposes, however, the size of each bin is 
exaggerated in Figure 2. 

Referring to Figures 1 and 2, at each of the ninety stops of the tomographic 
15 imaging device 30 along the 180' circumferential arc 31 about the object 28, 

gamma radiation emitted by the object 28 is received at the bins 42 of the 
detector 38. The collimator 36 serves to absorb any gamma radiation which is 
not perpendicular or normal to the plane of the detector 38. For example, a first 
gamma photon 44 shown in Figure 2 is not normal to the plane of the detector 
20 38, and is therefore absorbed by a wall of the collimator 36 without reaching any 

of the bins 42. A second gamma photon 46 travelling nomrial to the plane of the 
detector 38 passes through the collimator 36 unaffected, and is absorbed at one 
of the bins 42 of the detector. Other gamma photons 48 and 49 may miss the 
tomographic imaging device entirely, and are not detected. 

25 At each stop along the arc 31 , the tomographic imaging device produces signals 

representing a tomography scan image. More particularly, the signals represent 
data values, each data value indicating the number of gamma photons or 
"counts" received at a particular corresponding bin 42 at that respective stop. 
Thus, for each stop, a set of data values representing the numbers of counts 

30 received at respective bins is produced. These sets of data values for each stop 
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are transmitted from the tomographic imaging device 30 to the programmable 
device 26. It will be appreciated that as a result of the collimator, each 
tomography scan image produced at each of the stops of the tomographic 
imaging device represents a two-dimensional projection of the three- 
5 dimensional radioactivity in the object 28. 

Programmable device 

Referring to Figure 3, the programmable device is shown generally at 26. In this 
embodiment, the programmable device 26 is implemented in the desktop 
computer 34 shown in Figure 1. The programmable device includes a 

10 processor 50, which in this embodiment includes a microprocessor 52 in 

communication with a random access memory (RAM) 54, a storage memory 56 
and an input/output (I/O) port 58. In this specification, the tenns "processor" or 
"processor circuit" are used to include any programmable device or any circuit or 
combination of circuits capable of performing the functions described herein. 

15 Alternatively, therefore, the processor 50 need not be implemented in a desktop 

computer, and may include a combination of one or more microprocessors, 
microcontrollers, other integrated circuits, or logic gate arrays, for example. 
Other such variations will be appreciated by one of ordinary skill in the art upon 
reading this specification and are not considered to depart from the scope of the 

20 present invention. 

In this embodiment, the storage memory 56 includes a hard disk drive. 
Alternatively, however, it will be appreciated that any other type of memory may 
be substituted therefor. The storage memory 56 includes an input data store 60 
for storing sets of data values produced by the tomographic imaging device in 
25 respective two-dimensional input data matrices. In this embodiment, there will 

be ninety two-dimensional input data matrices 62, each input data matrix 
representing a respective tomography scan image produced by the tomography 
imaging device 30 at a particular one of the ninety stops or positions of the 
imaging device relative to the object. 
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The storage memory further includes a dynamic data store 64, for storing ninety 
three-dimensional dynamic data matrices 66, each dynamic data matrix 
representing the radioactivity of the voxels 40 of the object 28 at a particular 
time at which a corresponding tomography scan image was produced, or in 
other words, a three dimensional radioactivity image of the object 28 at a 
particular time during the tomography scan,. 

The storage memory 56 also includes a routines area 68 for storing a probability 
routine 70, a receive data routine 72, and a dynamic solution routine 74. 

The probability routine 70 directs the processor 50 to calculate probability 
coefficients representing probabilities that photons emitted from each voxel will 
arrive at a particular one of the bins 42 of the detector 38 at a particular time 
interval or detector position during the scan. 

The receive data routine 72 directs the processor 50 to receive signals 
representing data, or more particularly, representing successive tomography 
scan images of the object 28, and to store the data in the input data store 60 in 
the storage memory 56. 

The dynamic solution routine 74 directs the processor 50 to impose a shape 
constraint such as a linear inequality constraint to determine sets of dynamic 
data values from the input data, and to store such sets in matrices in the 
dynamic data store 64 in the storage memory 56. 

It will be appreciated, however, that the storage memory 56 is merely one 
example of a computer readable medium for providing computer readable 
instructions for directing a programmable device to implement the methods 
disclosed herein. Alternatively, the instructions or routines may be stored on a 
compact disc or floppy diskette, or on a ROM or EEPROM chip, for example, or 
on any other computer readable medium. Similarly, alternative means of 
producing a computer data signal embodied in a carrier wave for directing a 
programmable device to implement the methods disclosed herein would be 
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apparent to one of ordinary skill in the art upon reading this specification. Such 
variations are not considered to depart from the scope of the present invention. 

The storage memory 56 further includes a configuration data store 80 for storing 
configuration data, such as the probability coefficients calculated by the 
5 processor 50 under the direction of the probability routine 70, for example, for 

future use and/or re-use. 

The storage memory also includes a display store 81 for storing illumination 
matrices representing images of the object 28 at successive time intervals. 

Referring to Figures 1 and 3, the I/O port 58 is in communication with the data 
10 gathering device 22, and acts as a receiver for receiving the sets of values 

produced by the data gathering device representing measurements of the 
measurable property across the object 28 at respective measurement times, 
each set being associated with a respective measurement time. More 
particularly, in this embodiment, data representing successive tomography scan 
1 5 images across the object 28 detected at a SPECT imaging device at respective 

measurement times are received at the I/O port 58. The I/O port 58 is also in 
communication with the display device 35, and acts as a transmitter for 
transmitting a time varying graphical representation or values representing the 
measurable property to the display device for displaying a representation of the 
20 property. 

The RAM 54 includes a probability store 82, a calculation area 84, a counter 
register 85, a weighting factor store 86, a shape constraint store 87 and a 
display buffer 89. The probability store 82 is used to store a linear operator 
comprising probability coefficients calculated by the processor during execution 

25 of the probability routine 70. The calculation area 84 is used by the processor 

50 under the direction of the dynamic solution routine 74, to solve for the 
dynamic data matrices 66 as described in greater detail below. The counter 
register 85 is used to store a counter value representing the number of the 
tomographic scan currently in progress. The weighting factor store is used to 

30 store a weighting factor calculated by the processor during execution of the 
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dynamic solution routine 74, the weighting factor representing a variance of the 
input data stored in the input data store 60. The shape constraint store is used 
to store a shape constraint, such as a linear inequality constraint for example, to 
which the physical property of each voxel of the object is expected to confonn 
5 throughout the duration of the tomographic scan. The display buffer 89 is used 

to store data corresponding to the contents of the dynamic data matrices 66 
representing radioactivity images of the object 28 at successive times, for 
successive display of such images on the display device 35 shown in Figure 1. 

Probability Routine 

10 Referring to Figures 2 and 3, the probability routine 70 shown in Figure 3 

includes blocks of codes which direct the processor 50 to calculate probability 
coefficients to determine the values of a linear operator In this embodiment, 
hj and k denote the voxels 40, the bins 42 and the ninety stops along the arc 31 
respectively (1 <i< 27000, 1 <j < 3300 and 1 <k< 90), Each of the values Cp 

1 5 represents the probability that a gamma photon emitted from the voxel of the 

object 28 will arrive at the bin 42 of the detector 38 at time tjc, or in other words 
at the time when the tomographic imaging device 30 is at the stop of the 
ninety stops along the 180° circumferential arc 31 about the object 28. The 
probability coefficients Qjfc are based primarily on the geometry of the 

20 tomographic imaging device 30 and its physical arrangement and orientation 

relative to the object 28. However, the probability coefficients may also include 
probabilistic effects such as gamma ray scattering, attenuation and collimator 
blurring. 

The probability routine 70 may be or may include a conventional routine such as 
25 those which are known to those of ordinary skill in the art, and need not be 

described in further detail here. Alternatively, any routine capable of 
determining the probability coefficients may be substituted therefor. 



The probability routine further directs the processor 50 to store the coefficients 
of the linear operator Qjk in the probability store 82 of the RAM 54. Optionally, if 
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similar scans are to be performed again, the probability coefficients Qjk may also 
be stored in the configuration data store 80 in the storage memory 56, for future 
re-use. 

Receive Data Routine 

5 Referring to Figures 3 and 4, the receive data routine is shown generally at 72 in 

Figure 4, Generally, the receive data routine 72 includes blocks of codes which 
direct the processor 50 to receive a plurality of sets of values representing 
measurements of the property across the object 28 at respective measurement 
times, and to store the values in the storage memory 56. In this embodiment, 
10 the physical property of interest is radioactivity of the object 28 at respective 

times at which respective tomography scan images are produced by the 
tomographic imaging device 30, 

The receive data routine begins with a first block of codes 92 which directs the 
microprocessor 52 to set a counter A: = 1 in the counter register 85 in the RAM 
1 5 54, to correspond to the first time interval 4 or in other words, the first position or 

stop of the tomographic imaging device 30 along the circumferential arc 31 
during the tomography scan. 

Block 94 then directs the microprocessor 52 to receive data signals representing 
a tomographic image, and to store data values representing the tomographic 
20 image in an input data matrix 62 in the input data store 60 of the storage 

memory 56. 

Each input data matrix 62 has the same number of entries as the number of 
receptor bins 42 of the detector 38 of the tomographic imaging device. In this 
embodiment, therefore, each input data matrix 62 is a 50 row x 66 column 
25 matrix, with each entry corresponding to a particular respective one of the 50 x 

66 receptor bins 42. For numerical computational convenience, however, each 
of the 3300 receptor bins 42 is assigned a unique address number j where 1 <j 
< 3300, and each respective corresponding location in the input matrix 62 is 
assigned an identical address number. Thus, in the remainder of this 
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specification, when tine f bin or f location of an input matrix is referred to, it will 
be understood tliat the bins and corresponding input matrix locations each have 
addresses numbered 1-66 along a top row, 67-132 along a second row, and so 
on, down to a bottom row of bins and corresponding input matrix locations 
5 numbered 3235-3300. For example, the location (32,57) of the input matrix 62 

and the corresponding bin 42 located in the 32""* row, 57^ column of the array of 
receptor bins are referred to as the 2103'"^ 0=2103) input matrix location and the 
2103"^*^ bin 42 respectively. 

Thus, as data signals produced by the tomographic imaging device 30 are 
10 received at the I/O port 58, block 94 directs the microprocessor 52 to store data 

values representing the data signals in the input matrix 62, in a manner such 
that the data value stored in the location of the input matrix 62 represents the 
number of gamma ray counts detected at the corresponding receptor bin 42 
during the time interval. For example, the data value stored in the 2103'^ 
15 location (i.e. row 32, column 57) of the input matrix 62 represents the number 

of gamma ray counts detected at the corresponding 2103'"^ receptor bin 42 (i.e. 
the bin at the 32"^ row, 57^ column of the array of receptor bins 42 of the 
detector 38), during the time interval. 

Block 96 then directs the microprocessor 52 to determine whether the 
20 tomography scan is complete. In this embodiment, ninety separate tomographic 

scan images corresponding to ninety different positions of the tomographic 
imaging device 30 along the arc 31 are produced, resulting in ninety 
corresponding input data matrices 62 in the input data store 60 of the storage 
memory 56. If all ninety such input data matrices have not yet been generated, 
25 block 98 directs the microprocessor 52 to increment the counter k in the RAM 

54, and the microprocessor is directed back to block 94 to await receipt of the 
next set of signals representing the next tomographic scan image of the object. 
Processing continues in this manner until all ninety input data matrices have 
been generated, at which point the receive data routine is ended. 
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Dvnamic Solution Routine 
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Referring to Figures 3 and 5, the dynamic solution routine is sliown generally at 
74. Generally, the dynamic solution routine 74 includes blocks of codes which 
program the processor 50 to produce a plurality of sets of values representing 
the measurable property at a plurality of locations throughout the object at the 
5 respective measurement times, by minimizing a figure of merit function relating 

the values representing the measurements with the values representing the 
property at the plurality of locations, with a shape constraint imposed on the 
values representing the property at the plurality of locations. The dynamic 
solution routine achieves this by directing the processor to perform calculations 

10 by imposing an inequality constraint to determine dynamic data values from the 

data stored in the input data matrices. Each dynamic data value represents a 
physical property (radioactivity) of the object at a respective corresponding one 
of the object's voxels 40 at a respective corresponding time inten/al More 
particularly, the dynamic solution routine determines the radioactivity x of each 

1 5 voxel i at each time interval K denoted by Xi(t0 for k= (1 90). 

In principle, ignoring noise, the problem of solving ior Xi(tp) amounts to solving an 
equation of the form: 

i 

where Qjk is the linear operator comprising the probability coefficients, Xik 
20 denotes the dynamic data values Xi(t0 to be solved for and stored in the dynamic 

data matrices 66, and djk denotes the tomographic scan data stored in the input 
matrices 62. In practice, however, perfect solutions to equation (1) cannot be 
found, and reliance upon approximation techniques, such as least squares 
methods for example, is necessary. 

25 It will be appreciated that if the distribution of the radiopharmaceutical from voxel 

to voxel within the object were static, that is, constant over all ninety time 
intervals, then the radioactivity x of each voxel i would conform to the model 
^/(^a) = 4^~^' where X is the known decay constant of the radioisotope label. 
However, due to non-uniform flows of the radiopharmaceutical within the object, 
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each voxel i may have its own unique decay profile. Accordingly, some 
researchers have previously forced the radioactivity function of each voxel i to 
conform to more complex models, such as the following biexponential model, for 
example: 



where there are five unknowns for each voxel i, namely ^6 h Bu [Xi and Q. 

It will be further appreciated that a non-linear least squares (NLS) method is an 
example of a prior art approximafion method capable of deriving solutions of the 
form of equation (2). According to the NLS method, for each voxel i, the five 
parameters Au h Bi, and C, in equation (2) are solved by minimizing a figure- 
of-merit function ("merit function") of the following form: 




(3) 



subject to: 



Xi(t0 defined by the i basis functions of the form (2), i.e., 

;c,(?,) = 4e-"'"+5,e-''-'*+C, ; 



where: 



djk denotes the data value stored in Vnef location of the input 
data matrix 62 corresponding to the f receptor bin, representing 
the number of counts detected at receptor b'mj during time interval 



Qjk is the linear operator comprising the probability coefficients 
computed by the processor while executing the probability routine; 
and 
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(fju is a weighting factor, determined for example from tine 
variance of the receptor bin contents djk (generally based on a 
Poisson statistic). 

However, a merit function of the fomn of equation (3) subject to the basis 
functions (2) is not necessarily unimodal. Accordingly, it is possible for some 
minimization algorithms to get "stuck" in the vicinity of local minima, and 
therefore fail to find the best fit corresponding to the global minimum. In 
addition, the computations required to minimize a merit function of the form (3) 
subject to basis functions of the form (2) are resource-intensive and time- 
consuming, often taking as long as one to two hours to solve with contemporary 
desktop computers. Moreover, the biexponential model which is imposed on the 
basis functions (2) and therefore on the solutions may not necessarily apply in 
all situations. Extensions of this model including more exponential terms may 
suffer from similar difficulties. 

In contrast, according to the present invention, no assumption is made as to the 
precise form of the basis functions or the solutions Xi(t)^. Rather, all that is 
assumed is a shape constraint, which in this embodiment takes the form of a 
linear inequality constraint. For example, in some circumstances it may 
reasonably be expected that the activity at any given voxel 40 of the object 28 
will decrease over the course of the tomography scans, in such a case, the 
assumed shape constraint takes the form of the linear inequality constraint; 

X, it, ) > X, (t^ ) > X- ) . . . > X . (t^, ) > X. (t,, ) (4) 

Or, alternatively, in other situations only washin or uptake will occur during the 
course of the tomography scans, in which case it may be reasonable to assume 
that the shape constraint takes the fomn of the linear inequality constraint: 

X. (t, ) < X. {t^ ) < X, (?3 ) . . , < X, (t,, ) < X. (t,, ) (5 ) 

Other similar shape constraints may be readily envisioned to apply to certain 
situations. For example, in some cases it may be desirable to assume that Xiftf^) 
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increases then decreases representing washin followed by washout, in which 
case it might be assumed that the shape constraint takes the form of the linear 
inequality constraint: 

X, (t, )<x.(t,)'"< X. {t,^_, ) < X, {t,^ ) > X. (r^,, ) . . . > X, {t,, ) > X, ) (6) 
where position Pt may vary from voxel to voxel. 

In other cases it may be desirable to assume that Xi(tj^ decreases then 
increases. Similarly, it may be desirable to assume that Xi(t0 increases or 
decreases in either a concave or convex manner, without assuming any precise 
functional fomn of Xi(tk). Further illustrative examples of shape constraints are 
provided in (8) below. 

Thus, according to the present embodiment of the invention, a novel and 
inventive reconstruction method using Shape Constrained Least Squares 
(SCLS) analysis is provided. The solutions Xi(t;J are calculated by minimizing a 
merit function of the form (3) subject to a linear inequality constraint of the fomn 
(4), (5) or (6) for example, instead of minimizing the merit function subject to 
basis functions of the form (2). 

More particularly, in this embodiment the dynamic solution routine 74 directs the 
processor 50 to solve the following equation: 




(7) 



subject to: 



Xik^x(tO are numerical values satisfying a shape constraint, 



such as one of the following exemplary shape constraints: 



(8) 



(8a) x,^ > x,^ > . . . > > X,,, > 0 (decreasing) 



(8b) 0 < jc.j < < X.3 , . . < x.g, < (increasing) 
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(8c) 0 < x,^ <x,,^-< < x,^ > ^,^,1 . . . ^ X,,, > X,,, > 0 
(increasing to peak at Pf then decreasing) 

f8d)0<x. <x.--*<x and X. ... > x ,^ > x ,,^ > 0 
(increasing then decreasing with gap around peak Pi) 

(Be) 0 > X-, - X.2 > - X -3 > . . . > x-g, - and 
X.J >0,...x.9o >0 (convex) 

(8f) 0 < X,, - x,2 < x,2 - < . . . < x,g, - x,,o and 
x-i >0,...x-9o >0 (concave) 

(8g) 0 > x,j - x,2 > - x,3 > . . . > x,^^ ~ jc,^ and 

-^.-..^1 ^^^-.,-1 -^/f?,-2 -x,-9o (convex- 

concave with inflection point at 

wherein: 

R(x) is an optional regularizing term which may be omitted if 
desired, as discussed in greater detail below; and 

djk. Cijk and (/jk are defined as above in equation (3). 

Thus, rather than imposing any pre-determined functional fomn such as a 
biexponential form upon the solutions x/ifj^, the present embodiment of the 
invention presumes nothing about the form of Xift^ and simply solves for Xi(tk)= 
Xii ... Xi9o where each one of x/7 to Xfpo is simply a numerical value representing 
the radioactivity at voxel i at measurement times ti to % respectively. 

It will be appreciated that minimizing the merit function (7) subject to a shape 
constraint such as those listed in (8) for example, presents an optimization 
problem where the inequality constraints are linear in the unknown variables Xik 
= x(tj^, representing the measurable property for each of the locations / 
throughout the object at the times k. Effectively, therefore, in this embodiment 
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equation (8) defines 27,000 respective linear basis functions Xi(tO= xn . . . Xi9o , 
each linear basis function fiaving ninety degrees of freedom. 

Referring to Figures 3 and 5, tiie dynamic solution routine is shown generally at 
74 in Figure 5. The dynamic solution routine 74 includes a plurality of blocks of 
5 codes including a first block of codes 100, which directs the processor 50 to 

calculate the weighting factor a/^ refen-ed to above in equations (3) and (7). 
Essentially, the weighting factor a/^ represents a covariance matrix and is used 
to account for the expected measurement noise. Block 100 first directs the 
microprocessor 52 to read the contents djk of the ninety input data matrices 62, 

10 each input data matrix containing 3300 values representing the number of 

gamma ray counts received at each / one of the 3300 bins 42 during the k*^ one 
of the ninety time intervals. Block 100 then directs the microprocessor to 
calculate the values of the weighting factor a/^ by determining a variance of the 
contents djk of the input data matrices. More particularly, in this embodiment the 

1 5 values of the weighting factor a/^ are determined from the Poisson statistics of 

the input data values djk. For example, the values of the weighting factor ajk'^ 
may be calculated such that ajk = max(djk, e) for a threshold value e>0 to avoid 
exceedingly small values o/ if the data djk are approximately Poisson 
distributed. Other methods for producing the weighting factor would be 

20 apparent to one of ordinary skill in the art upon reading this specification, and 

such alternative methods are not considered to depart from the scope of the 
invention. 

Block 100 further directs the microprocessor to store the weighting factor ay/ in 
the weighting factor store 86 of the F?AM 54. Optionally, if further analysis is to 
25 be repeated for the same input data, block 100 may further direct the processor 

to store the weighting factor in the configuration data area 80 of the storage 
memory 56. 

The microprocessor 52 is then directed to block 102, which directs it to 
detennine an appropriate shape constraint for the i solutions Xift/;). In this 
30 embodiment, block 102 prompts a user of the programmable device 26 to enter 
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the applicable shape constraint at the user input device 33 shown in Figure 1. 
For example, in response to receiving user input of "DECREASING", the 
microprocessor determines that the linear inequality constraint (4) applies, and 
stores the constraint (4) in the shape constraint store 87 of the RAM 54. 
5 Alternatively, the microprocessor may be directed to read a pre-determined 

shape constraint previously stored in the configuration data area 80 of the 
storage memory and copy this stored constraint to the shape constraint store 
87. 

Block 104 then directs the microprocessor 52 to minimize the merit function (7) 
10 above, subject to the linear basis functions comprising shape constraints of the 

form (8) for example. To do this, the microprocessor 52 is directed to execute a 
math program, or more particularly a minimization subroutine (not shown) on the 
values representing the measurements, to minimize the sum of squares (7) of 
the difference between the product HiCykXik of the linear operator and the values 
15 representing the property, and the values djk representing the measurements of 

the property. 

In this embodiment, the minimization subroutine may include a conventional 
Constrained Optimization subroutine. However, block 104 directs the 
microprocessor to execute the subroutine in an unconventional manner, by 
20 substituting linear basis functions comprising a shape constraint of the form (8) 

stored in the shape constraint store 87, for conventional basis functions which 
are normally used in such a minimization subroutine. 

More particularly, in this embodiment a conventional bound constraint routine is 
modified to minimize equation (7) subject to shape constraints, such as linear 

25 inequality constraints. Such conventional bound constraint software may 

include an L-BFGS routine described in the publication by D. Liu and J. 
Nocedal, "On the Limited Memory BFGS Method for Large Scale Optimization", 
Math. Prog. B. 1989 at pp. 503-528, for example, or alternatively may include a 
GPCG routine by Rolf Felkel of the Numerical Analysis Group, Mathematical 

30 Department, Technical University Darmstadt, for example. 
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Generally, such conventional bound constraint software has been previously 
used to minimize nonlinear objective functions subject to bound constraints on 
the unknown variables Xi, that is, h < xt < Ui for given boundary values 

According to the present embodiment of the invention, block 104 directs the 
5 microprocessor to perform a change of variables to effectively use such 

software to achieve minimization of equation (7) subject to shape constraints. 
For example, where the shape constraint stored in RAM is of the form (8a), (8b), 
(8c) or (8d), block 104 directs the microprocessor to perfonn the following 
change of variables: 

10 set (9) 

[■^190 ~ '"■190 

Similarly, where the shape constraint stored in the RAM is of the fonn (8e), (8f) 
or (8g), block 104 directs the microprocessor to perform the following change of 
variables: 

set xl = 4 - x;^, =x^- 2x„^, + (1 0) 

1 5 It Will be appreciated that the changes of variables (9) or (1 0) serve to transform 

the shape constraints on xik of the form (8a) - (8d) or (8e) - (8f) into bound 
constraints on x'ik or respectively. Altematively, other transformations of 
other types of shape constraints will be apparent to those of ordinary skill in the 
art upon reading this specification and are not considered to depart from the 

20 scope of the invention. 

The shape constraint previously stored in RAM at block 102 is effectively pre- 
specified to the math program or minimization routine invoked at block 104. 
Thus, block 104 directs the microprocessor to read the pre-specified shape 
constraint, transfonn the shape constraint into a bound constraint according to 
25 (9) or (10) or other transformations, to execute the conventional bound 

constraint minimization routine subject to the change in variables, then to 
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transform back to Xik, to effectively minimize equation (7) subject to the pre- 
specified sliape constraint. 

Optionally, block 104 may direct the processor to stabilize the minimization 
problem by adding various regularizing terms to the quadratic objective in 
equation (7). This would entail minimizing equation (7) wherein at least some of 
the regularizing terms R(x) are non-zero. Examples of preferable regularization 
terms include the following: 

ik 

wherein R(x) is a Tychonoff type (weighted spatial gradient) 
regularizer, and wherein f, ^ and f are voxels adjacent to voxel 1, 
such as an upper neighbour, a right hand neighbour and a forward 
neighbour, for example; 

W = Z ^iki^ik.^-^ikf (11b) 
ik 

wherein R(x) is a Tychonoff type (weighted gradient in time) 
regularizer; 

k i 

wherein R(x) is a Fourier domain type (spatial Fourier high pass 
filter) regularizer, using the spatial Fourier transform of each of the 
image frames; the weights are close to zero for low spatial 
frequencies / and close to 1 for high spatial frequencies i ; 

/ k 

wherein R(x) is a Fourier domain type (temporal Fourier high pass 
filter) regularizer, using the temporal 1 D Fourier transfonn for each 
of the time series xa ... xm belonging to voxel i\ the weights are 
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close to zero for low temporal frequencies k and close to 1 for high 
temporal frequencies A: ; or 



(11e) 



wherein R(x) is a Physically motivated (constant flow rate) 
5 regularizer which gives preference to time profiles xa ... Xi9o 

exhibiting a constant flow rate of the form Ae^\ and therefore 
penalizes irregular flow profiles; this regularizer acts in the 
temporal domain on each pixel separately. 

Although regularizing terms may be omitted entirely if desired, it has been found 
10 that streaking artifacts, which may otherwise arise in images corresponding to 

early reconstruction times, tend to subside upon incorporation of such 
regularization. The Tychonoff and Fourier regularization methods have been 
found to provide slightly better numerical performance than the physically 
motivated method. 

15 The microprocessor thus combines minimization and regularization to produce, 

for each voxel / of the voxels 40 of the object 28, a solution of the form (8) 
described above. 

Block 104 further directs the microprocessor 52 to store the solutions Xi(t0 in 
corresponding locations in the dynamic data store 64 in the storage memory 56. 

20 As noted above, in this embodiment the dynamic data store 64 contains ninety 

dynamic data matrices 66, each matrix corresponding to a particular time 
interval k at which a corresponding tomographic image was produced. More 
particularly, in this embodiment, each of the dynamic data matrices 66 in the 
dynamic data store is a 30 x 30 x 30 three-dimensional matrix for storing 27,000 

25 dynamic data values, each of the data values corresponding to a particular one 

of the 27,000 voxels 40 of the object 28 shown in Figure 2. Thus, for each 
voxel z, the microprocessor produces ninety dynamic data values Xik for times 
tk=(ti.,.t9o). The microprocessor first stores xa in the location of the first 
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dynamic data matrix, then stores Xi2 in the location of the second dynamic 
data matrix, and so on, until all ninety of the values Xik corresponding to the 
voxel have been stored in respective dynamic data matrices. The 
microprocessor continues in this manner until ninety such values have been 
5 stored in the dynamic data matrices representing the number of counts at each 

of the 27,000 voxels i at each of the ninety time intervals L 

It will be appreciated that because the present invention need not be 
implemented in a desktop computer, the particular format of the signals 
produced by the processor 50 at block 104 is unimportant. Generally, any 
1 0 alternative type of signal representing the plurality of sets of values representing 

the property at the plurality of locations throughout the object is not considered 
to depart from the scope of the present invention, provided the signal is 
produced according to the inventive methods disclosed herein or equivalents 
thereof, 

1 5 The processor is then directed to block 106, which directs the microprocessor to 

produce a time-varying graphical representation of the object 28, to represent a 
change of the measurable property over time. The processor is then directed to 
use the representation to control the display device to display the 
representation, or more particularly, to produce successive images representing 

20 the property at successive instants in time, in response to the representation of 

the property. This is effectively achieved by successively displaying 
representations of the respective dynamic data matrices 66. 

It will be appreciated that a number of display options are available. In the 
present embodiment, the microprocessor is directed to produce a two- or three- 

25 dimensional projectional representation of the three-dimensional object 28 and 

to display the projectional representation on the display device 35. Block 106 
directs the microprocessor to execute a projectional display subroutine (not 
shown) to produce A: two-dimensional illumination matrices, such that each value 
in each illumination matrix represents an intensity of illumination of a 

30 corresponding pixel on the display. 
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Each of the k illumination matrices represents a view at a successive time frame 
of the object 28, possibly from a different successive spatial point of view, so 
that a successive display of the contents of the k illumination matrices is akin to 
viewing a short movie of the 3D object in which the camera slowly moves from a 
5 first point of view to a second point of view. To achieve this, block 106 directs 

the microprocessor to produce each of the k illumination matrices to represent a 
projection of the corresponding three-dimensional dynamic data matrix onto a 
hypothetical two-dimensional an-ay, with the hypothetical two-dimensional array 
incrementally moving along a radial arc about the three-dimensional dynamic 
1 0 data matrix, analogously to the motion of the tomographic imaging device about 

the object, so that each successive illumination matrix represents a view of the 
object from an incrementally different point of view than the previous view. The 
"views" represent snapshots of the dynamic object at successive time 
increments. 

15 For example, referring to Figure 6, a schematic representation of a dynamic 

data matrix is shown generally at 66, Each of the locations in the dynamic data 
matrix has an identifying number f, 1 <i< 27,000, uniquely corresponding to a 
particular one of the 27,000 voxels 40 of the object 28 shown in Figure 2. It will 
be appreciated, for example, that in one possible view at time k=l, voxel 

20 numbers 1 to 30 of the 30^ voxels 40 might be along a line normal to the 

hypothetical two-dimensional array, and accordingly, their values would be 
projected into the same location of the 1®* illumination matrix, so that a single 
illumination value in the 1®* illumination matrix corresponding to a single pixel 
would represent the summed radioactivity of voxels 1 to 30, or in other words 

30 

25 ^^ix ■ However, in another possible view at time k=46, for example, voxel 

numbers 1, 901, 1801, 2701, ... 26101 might now be along the line normal to 
the hypothetical two-dimensional array, in which case a single illumination value 
in the 46*^ illumination matrix corresponding to a single display pixel would 
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contain the summed value '^x^i^goo^,-_i^y^46 representing the projection of the 
radioactivity of voxels along the normal line. 

Referring back to Figures 3 and 5, block 106 further directs the microprocessor 
52 to normalize the illumination matrices, if necessary to adjust for the range of 
5 pixel illumination available on the display device 35, and to store the illumination 

matrices in the display buffer 89 in the RAM 54. In addition, the microprcicessor 
is directed to store the illumination matrices in the display store 81 of the storage 
memory 56. Once all of the illumination matrices have been produced and 
stored in this manner, block 106 directs the microprocessor to successively 
10 display the illumination matrices, to effectively produce a visual representation of 

the object in response to the plurality of sets of values representing radioactivity. 
The visual representation includes a movie representing the radioactivity of the 
voxels 40 of the object 28 over the ninety time intervals. 

The foregoing projectional display method is merely one example of a way in 

15 which the processor may be programmed to produce the successive images to 

depict a three dimensional representation of the property in the object which 
varies according to changes in the property over time. Alternatively, other ways 
of producing a visual representation of the object in response to a time varying 
graphical representation produced by the processor may be apparent to one of 

20 ordinary skill in the art upon reading this specification. For example, different 

display colours may represent different intensities of radioactivity at 
corresponding voxels of the object. Or, as a further example, radioactive local 
regions may be displayed in a manner such that the size of the local region is 
exaggerated proportionally to the radioactivity of the region. Alternatively, the 

25 contents of the dynamic data matrices may be used to generate true three- 

dimensional images comprising left and right perspective views, preferably from 
a gradually shifting point of view, and the display device 35 shown in Figure 1 
may be modified to provide for synchronized viewing glasses to enable three- 
dimensional viewing, for example. A 2-D slice of the object may be displayed, or 

30 a 3D dynamic movie may be displayed by showing iso-surfaces of the 3D object 
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at each time, possibly from spatially shifting points of view. Or, volume 
rendering may be used to display transparent 3D objects. Alternatively, it may 
be desirable to observe the time profiles for one or more individual voxels, 
displayed as static graphical curves, for example. More generally, the 
representation produced at block 106 is not limited to a display and may include 
other fomns of representations. For example, the representation may include 
any fomn of kinetic modelling. Other similar variations will be apparent to one of 
ordinary skill in the art upon reading this specification, and are not considered to 
depart from the scope of the invention. 

Generally, it has been found that the foregoing reconstruction method of Shape 
Constrained Least Squares (SCLS) analysis provides numerous advantages 
over the prior art. Generally, the singular value distribution of the matrices Ojk^ 
djk encountered in many practical situations allows for the constrained 
optimization subroutine to achieve fast convergence. Thus, in contrast with 
conventional Nonlinear Least Squares methods or the Factor Analysis of 
Dynamic Stnjctures method which often require one to two hours or more of 
processing time, embodiments of the present invention often solve for the 
dynamic data Xi(t]^ for all voxels i over all time intervals k within ten minutes, for 
example. 

In addition, since no assumption is made as to the precise functional form of the 
solutions Xi(t0, the present invention does not sacrifice experimental accuracy by 
attempting to force the data to conform to a biexponential or other standard 
model. 

For each voxel i, the inventive SCLS reconstruction method described herein 
effectively allows for an activity profile with a much larger number of degrees of 
freedom, such as ninety degrees of freedom where the object is scanned at 
ninety successive time intervals, for example. Although conventional wisdom 
suggests that increasing the degrees of freedom generally increases the risk of 
producing an under-fitted model, experiments have generally confirmed that the 
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SCLS method disclosed herein works well and amounts to a significant 
improvement over existing NLS, FADS or fast rotation methods. 

Alternatives : 

Other Minimization Routines: 

5 Embodiments of the present invention may employ constrained minimization 

routines other than constrained least squares. 

For example, referring back to Figures 3 and 5, in a second embodiment of the 
invention the dynamic solution routine 74 shown in Figure 3 may include a 
Maximum Likelihood (ML) minimization subroutine. Block 104 in Figure 5 may 
1 0 be modified to direct the processor 50 to solve an equation of the following form: 



Minimize \(I)(^) = Y.Y. S ^iA^^ik - dj, log S Q^^- 



3300 90 f 27000 f 27000 




ik 




(12) 



subject to: 



Xik = x^t0 being defined as a linear basis function of the shape- 
constrained fonm (8). 
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Block 104 may further be modified to combine the maximum likelihood 
minimization subroutine with regularization, as previously described. 
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Similarly, it will be appreciated that according to a third embodiment of the 
invention, the dynamic solution routine 74 shown in Figure 3 may include an 
Expectation Maximization (EM) subroutine, which may include an OS-EM or 
Bayesian EM routine, for example. Block 104 in Figure 5 may be modified to 
direct the processor 50 to execute the EM subroutine, preferably with 
regularization as described above, to solve equations analogous to equation (7) 
subject to (8) above, which may be easily derived by one of ordinary skill in the 
art upon reading this specification. 
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According to a fourth embodiment of the invention, prior information about the 
unknown tracer dynamics in the body may be used, if available. As most 
physiological models may be represented by compartmental models, it may be 
desirable to include such a model in a given situation. For example, if the 
dynamic to be visualized is known to be accurately represented by a two 
compartment model, a biexponential function of the form (2) may be correct. 
Rather than imposing the biexponential functional form, the present embodiment 
of the invention seeks reconstructions which are close to satisfying a difference 
equation of the following form; 

It will be appreciated that a set of values x^k = Xi(t^ satisfying equation (2) will 
also satisfy a difference equation of the form (14). 

Thus, referring back to Figures 3 and 5, in this embodiment, block 104 of the 
dynamic solution routine 74 is modified to direct the processor 50 to solve the 
following equation: 

Minimize: 



j.k \ i J i k 



(15) 



subject to: 

Xik=Xt{t0 satisfying shape constraints of the form (8); 

where R(x) comprises regularizing terms, such as those described 
above, for example; and 

^ is a penalty parameter. 



The processor 50 is directed to minimize the merit function F(x,a,p,Y) with 
respect to the variables (x,a,(i,y) simultaneously. 
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Modified block 104 first directs the processor to begin by setting an iteration 
value n=0 and a first penalty parameter K^O in the calculation area of the RAM 
54, and to obtain an estimate forx/> The processor is then directed to use the 
estimate for Xi^k to estimate the three parameters au A and y, for each voxel i. 
5 Modified block 104 then directs the processor to select a new incremented 

iteration value n>0, and a new penalty parameter K and to repeat the above 
estimation procedure using the previous estimates to obtain an improved 
estimate oi Xi^h Block 104 directs the processor to continue incrementing the 
iteration value n and adjusting the penalty parameter K to obtain successively 

10 improving estimates of x,;^ In practice, only a few iterations, usually of the order 

of five to ten iterations, are required for the successive estimates of to 
converge. Thus, in this embodiment, block 104 directs the processor to perform 
ten iterations. Generally, the choice of each successive penalty parameter K for 
each iteration is heuristic, and must be "tuned" to the particular application, by 

15 trial and error. Successive penalty parameters are non-decreasing, and are 

subject to the constraints that if the penalty parameter is too large the iteration 
will fail, but if it is too small the penalty parameter will have no influence on the 
iteration. In this embodiment, block 104 directs the processor to select the ten 
successive penalty parameters as ^ = {0, 2, 5, 10, 15, 20, 20, 20, 20 and 20} for 

20 the ten iterations respectively. However, other suitable choices may be 

apparent to those of ordinary skill in the art upon reading this specification. 

Alternatively, rather than perfomning a fixed number of iterations, block 104 may 
be modified to compare a given estimate for x^k to a previous estimate and 
cease iterating when the difference falls below a pre-defined threshold, 
25 indicating substantial convergence. It has been found that the converged 

estimates for Xi^k generally agree well with biexponential results produced by 
solving equation (3) subject to (2). 

Although the measurable property included radioactivity in the foregoing 
embodiments, it will be appreciated that the embodiments of the present 
30 invention may be applied to other measurable properties in other contexts. 
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While specific embodiments of tiie invention have been described and 
Illustrated, such embodiments should be considered illustrative of the invention 
only and not as limiting the invention as construed in accordance with the 
accompanying claims. 



